% estimate impulse responses to government spending shock 
% this version: may 2011

clear all; clc; close all;

load us_quarterly       % generated by getdata.m


saveirf      = 1;

rand('state',1);
nptsvar     = 30;       % number of points of IRF in VAR
teststatio  = 0;        % tests wether VAR is explosive
confidence  = 1.645;    % (1.645: 90 percent); % or 1.96 %(95 percent);
ndraws      = 100;      % number of replications in bootstrap
nlagq       = 4;        % number of lags
beg_samp    = (82-47)*4+1; 
end_samp    = (107-47)*4+4;

disp(['Sample: ' num2str(sampleq(beg_samp+nlagq)) ' - '  num2str(sampleq(end_samp)) ' (dependent variables)  '  ])
gy = mean(g(beg_samp:end_samp)./y(beg_samp:end_samp));
cy = mean(c(beg_samp:end_samp)./y(beg_samp:end_samp));
iy = mean(inv(beg_samp:end_samp)./y(beg_samp:end_samp));

    
for sensi =1:4;

    if sensi == 4
        beg_samp    = (74-47)*4+1; 
    end
    for jrun = 1:3;

        if jrun==1;
            namshort = char('news','g','y','c','lr','rx','dy','inf');       
            xdata = [news log(g) log(y) log(c) lrr log(rx) dy  inflation];   
            if sensi==3 || sensi == 4;
                xdata = [news log(g) log(y) log(c) tbill log(rx) dy  inflation];   
            end
            yshare = [1 gy 1 1 1 1 1 1];
        elseif jrun ==2;
            namshort = char('news','g','y','i','lr','rx','dy','inf');       
            xdata = [news log(g) log(y) log(inv)  lrr log(rx) dy inflation];
            if sensi==3 || sensi == 4;
                xdata = [news log(g) log(y) log(inv)  tbill log(rx) dy inflation];
            end
            yshare = [1 gy 1 iy 1 1 1 1];
        elseif jrun ==3;
            namshort = char('news','g','y','nx','lr','rx','dy','inf');      
            xdata = [news log(g) log(y) nx  lrr log(rx) dy inflation];
            if sensi==3 || sensi == 4;
                xdata = [news log(g) log(y) nx tbill log(rx) dy inflation];
            end
            yshare = [1 gy 1 cy 1 1 1 1];
        end

    
    %% specify deterministics
    [obs nvar]=size(xdata);
    exo = []; 
    for ik=1:obs; ltrend(ik)=ik; end
    for ik=1:obs; exosq(ik)=ik^2; end
    exo = [ones(obs,1) ltrend' ];
    if sensi == 2
        exo = [ones(obs,1) ltrend' exosq'];
    end
   
    
    exo = exo(beg_samp:end_samp,:);
    samplelength=length(xdata(beg_samp:end_samp,:));
    
    
    %% estimate impulse responses
    [impzout,erzout,azeroout,a0betazout,var_explode]=...
        estimateIRF(xdata(beg_samp:end_samp,:),nlagq,nptsvar,teststatio,exo);       

    irf1(:,:,sensi,jrun) =  impzout(:,:,1);
   
    %% bootstrap IRF
        simdout=simdata(xdata(beg_samp:end_samp,:),...
        nlagq,nptsvar,ndraws,impzout,erzout,azeroout,a0betazout,exo);
        simimpzout=[];
        for zx=1:ndraws;
            [simimpzout,erzout,azeroout,a0betazout,var_explode] = estimateIRF(simdout(:,:,zx),...
            nlagq,nptsvar,teststatio,exo);                   
            irf1_mc(zx,:,:) = squeeze(simimpzout(:,:,1));
        end
    
        irfse1(:,:,sensi,jrun) = squeeze(std(irf1_mc(:,:,:),0,1));

    end
    
    
end

for jrun=1:3;
     if jrun==1;
                namshort = char('news','g','y','c','lr','rx','dy','inf');       
                yshare = [1 gy 1 cy 1 1 1 1];
            elseif jrun ==2;
                namshort = char('news','g','y','i','lr','rx','dy','inf');       
                yshare = [1 gy 1 iy 1 1 1 1];
            elseif jrun ==3;
                namshort = char('news','g','y','nx','lr','rx','dy','inf');       
                yshare = [1 gy 1 1 1 1 1 1];
     end    
        %% plot
    time = [0:nptsvar-1];
    jj=1;
    figure
    for zx=1:nvar;
        subplot(3,3,jj)
        mycolors = [1 1 1; .85 .85 .85];            
        colormap(mycolors);
         
        rescale = max(irf1(1:10,2,1,jrun))*(yshare(2)/yshare(zx));
        point = irf1(:,zx,1,jrun)./rescale;
        sepoint = irfse1(:,zx,1,jrun)./rescale;
        ub = point+confidence*sepoint;
        lb = point-confidence*sepoint;  
        hold on;   
        gx=[lb ub-lb];
        area(time,gx,min(lb),'EdgeColor','none');
        plot(time,zeros(nptsvar,1),'k--','LineWidth',1);
        %plot(time,point,'b-','LineWidth',1);
        
        sensi = 2;
        rescale = max(irf1(1:10,2,sensi,jrun))*(yshare(2)/yshare(zx));
        point = irf1(:,zx,sensi,jrun)./rescale;
        plot(time,point,'r--','LineWidth',1);
        
          sensi = 3;
            rescale = max(irf1(1:10,2,sensi,jrun))*(yshare(2)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'k-.','LineWidth',1);

            sensi = 4;
            rescale = max(irf1(1:10,2,sensi,jrun))*(yshare(2)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'b-','LineWidth',1,'MarkerSize',2);

            
            axis([-.2 nptsvar-.8  min([1.5*min(lb) 0]) 1.5*max(ub)]);box on;       
        
        title(deblank(namshort(zx,:)));
        jj=jj+1;
    end
end
    

  if saveirf
        
      
      
        disp('writing figures to files...')
      
        for jrun=1:3;
     if jrun==1;
                namshort = char('news','g','y','c','lr','rx','dy','inf');       
                yshare = [1 gy 1 cy 1 1 1 1];
            elseif jrun ==2;
                namshort = char('news','g','y','i','lr','rx','dy','inf');       
                yshare = [1 gy 1 iy 1 1 1 1];
            elseif jrun ==3;
                namshort = char('news','g','y','nx','lr','rx','dy','inf');       
                yshare = [1 gy 1 1 1 1 1 1];
     end    
        
        for zx=1:nvar;            
            figure
            set(gcf,'Visible','off')
            mycolors = [1 1 1; .85 .85 .85];            
            colormap(mycolors);
            rescale = max(irf1(1:10,2,1,jrun))*(yshare(2)/yshare(zx));
            point = irf1(:,zx,1,jrun)./rescale;
            sepoint = irfse1(:,zx,1,jrun)./rescale;
            ub = point+confidence*sepoint;
            lb = point-confidence*sepoint;  

            hold on;   
            gx=[lb ub-lb];
            area(time,gx,min(lb),'EdgeColor','none');
            plot(time,zeros(nptsvar,1),'k--','LineWidth',1);
            %plot(time,point,'b-','LineWidth',1);
            axis([-.2 nptsvar-.8  min([1.5*min(lb) 0]) 1.5*max(ub)]);box on;       
            set(gca,'FontSize',25)
            box on;
            
               sensi = 2;
        rescale = max(irf1(1:10,2,sensi,jrun))*(yshare(2)/yshare(zx));
        point = irf1(:,zx,sensi,jrun)./rescale;
        plot(time,point,'r--','LineWidth',4);
        
          sensi = 3;
            rescale = max(irf1(1:10,2,sensi,jrun))*(yshare(2)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'k-.','LineWidth',4);

            sensi = 4;
            rescale = max(irf1(1:10,2,sensi,jrun))*(yshare(2)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'b-','LineWidth',4,'MarkerSize',2);
            
            saveas(gcf,['..\figures\Ram_sensi_' num2str(jrun) '_' deblank(namshort(zx,:)) '.eps'],'psc2')
        end

        end
  end

